function dy = demean_group(y, gp, n_gp, ind_gp, iS)

cumsum_iS = cumsum(y(iS));

y_mean = diff([0; cumsum_iS(ind_gp)])./n_gp;

dy = y - y_mean(gp, :);